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I . INTRODUCTION 


There are in existence several user-oriented digital computer pro- 
grams for the static analysis of shells of revolution. A detailed 
discussion of most of these programs is given in Ref. 1. Of particular 
interest here is the program developed by Ball [Ref. 2] for the geo- 
metrically nonlinear analysis of arbitrarily loaded shells of revolution 
This program is an equilibrium program; that is, it solves for the dis- 
placement and stress resultant fields for an arbitrary loading condition 
Since geometric nonlinearities are included, the magnitude of load that 
leads to a condition of instability can be determined. 

The utility of the program would be considerably enhanced if it 
could be used to determine bifurcation buckling loads and the behavior 
of the shell in the vicinity of the bifurcation load. This latter 
feature is often referred to as the imperfection sensitivity of the 
shell to the load. As a consequence, the objective of this study was 
to use the computer program to examine the buckling behavior of several 
shells subjected to axisymmetric and nearly axisymmetric loads. It was 
anticipated that an examination of the effects of the small asymmetric 
perturbations upon the stability of the shell would disclose the bi- 
furcation buckling load and provide a quantitative evaluation of the 
imperfection sensitivity of the shell to the load. 


II . DESCRIPTION OF THE COMPUTER PROGRAM 


The computer program used in this study is a modified version of the 
program described in Ref. 2. The program is capable of analyzing thin, 
linearly elastic shells of revolution subjected to arbitrary loads and 
temperature distributions. It solves Sanders' nonlinear field equations 
for the conditions of small strains and moderately small rotations. 

A. COORDINATE SYSTEM 

The coordinate system used to locate points within the shell is s, 
e and Q, where s is the meridional distance along a reference surface, 

Q is the normal to the surface and e is the circumferential angle 
measured about the axis of revolution. 

B. METHOD OF SOLUTION 

The governing partial differential equations are reduced to four 
second order ordinary differential equations by expanding all dependent 
variables in trigonometric series in the e direction. Each of the non- 
linear terms, which are products of series, is further expanded into a 
single series. The coefficients of like trigonometric arguments are 
then grouped to form coupled sets of ordinary differential equations. 

The sets of differential equations are uncoupled by treating the non- 
linear terms as known quantities or pseudo loads. A finite difference 
formulation is employed for the meridional derivatives of the variables. 
This leads to sets of algebraic equations of the form 

[ K Cn)] {,(»)}. { q (n,} + { q („)J m 
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The square matrix 


K°] * 


is a banded matrix associated with the linear 


portion of the governing equations . The unknown variables in the column 


matrix -iZ 


ix £ 


.(n)-l 


th 

are the coefficients of the n term in the series expan- 


sions of U, V, W and M at each finite difference station. The variables 

s 

U, V and W are the displacements in the s, e and Q directions respectively, 
and M is the meridional bending moment per unit length. The applied loads 


in the n*” mode are contained in the column matrix ■jq' n ^|, and the co 
matrix {q^} 


lumn 


contains all of the nonlinear terms. There is one such 
matrix equation for every value of n used in the trigonometric expansions 
of the dependent variables. 

An elimination procedure is used to solve each set of equations for 
the actual loads plus an estimate of the pseudo loads. The nonlinear 
terms are then recalculated, using the new solution, and entered as 
revised pseudo loads. All of the sets of equations are then solved 
again, and this procedure is repeated until two consecutive solutions in 
each mode differ by less than a specified amount. This method of solution 
is referred to as the Jacobi method. 

A load-displacement history is obtained by incrementing the applied 
loads in equal load steps until the number of iterations required to 
achieve convergence in any mode exceeds a prescribed limit . The load 
increment is then reduced by a factor of five, and the procedure con- 
tinues until a prescribed number of load steps or a prescribed number 
of load step reductions have been made. 

Since the method of solution is based on a nonlinear pseudo load 
approach, the shell reacts equally, in a linear fashion, to any change 
in either the applied load or the pseudo load. Thus, failure of the 
solution to converge in any mode can be attributed to two types of non- 
linear behavior. Both types are illustrated in Fig. 1. The existence 
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of a maximum or an inflection point on the softening load-deflection 
curve A represents a type of behavior for which a solution can be 
obtained only below the point of zero or nearly zero slope. On the 
other hand, the existence of a stiffening nonlinearity, as illustrated 
by curve B of Fig. 1, can also cause a convergence failure whenever the 
slope becomes too steep. Thus, in general, it is necessary to examine 
the load-displacement behavior of the shell in order to determine the 
cause of the convergence failure. 

C . MODIFIED PROGRAM 

The program used in this study was a modified version containing 
revisions made by personnel at the NASA Langley Research Center and by 
the original author. The main difference between the modified program 
and the version given in Ref. 2 is the manner in which core storage is 
allocated for the solution vector. The solution vector was changed from 
a three dimensional array to a two dimensional array, allowing any com- 
bination of meridional and circumferential unknowns to be specified. 

The product of the number of meridional stations and the number of terms 
in the trigonometric expansions must be less than 201, and the maximum 
number of terms in the trigonometric expansions must be less than 11. 

The difference in storage allocation necessitated changes in the 
method of handling the geometry, the inplane and bending stiffnesses, 
and the loads. The subroutine STIF was replaced by the subroutine BDB ; 
subroutine PTLOAD was replaced by two subroutines, PLOAD and TLOAD ; and 
subroutine ACOEF was eliminated. In addition, two subroutines, HJ and 
EFG, given in Ref. 3, were substituted for the original subroutines HJ 
and EFG. 
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LOAD 



Fig. 1 Types of nonlinear behavior that 
have an effect upon convergence. 
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An additional change has been made in the test for convergence. 
Convergence was set such that the difference in two successive modal 
solutions, at every station, must be less than the convergence criterion 
times the maximum solution in that mode, provided the maximum solution 
is greater than 10 ^ . The test for convergence is not required when the 
maximum solution is less than 10 ^ . This change was made to provide 
converged, accurate solutions in fewer iterations. 

Finally, the subroutine OUTPUT has been modified in order to present 
the data in more compact form; the COMMON and DIMENSION statements have 
been changed to allow the compilation of the program in any order; and 
several "bugs" were detected and eliminated. 

The modified program occupies approximately 146,000 bytes or about 
36,500 words on the IBM 360/67. 

D. PROGRAM INPUTS 

The user-prepared subroutines BDB, GEOM and PLOAD, describing the 
shell stiffness, geometry and applied loads, for the shells covered in 
this study are given in Appendix A. The important parameters that can 
be varied for a given shell geometry are; number of points, number of 
modes required in the solution, number of modes used to describe the 
load, boundary conditions, loading condition, convergence criterion, 
load step size, number of load step reductions and number of iterations. 
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III. APPLICATIONS 


The program was used to solve for the response of several shells 
subjected to axisymmetric and nearly axisymmetric loads. Figure 2 
presents three load-displacement curves that are typical of the types 
of axisymmetric behavior predicted by the program for a shell that is 
imperfection sensitive. 

A. SNAP BUCKLING 

Axisymmetric snap buckling is illustrated by curve E of Fig. 2, and 
the application of the program to determine the snap buckling load is 
fairly straight forward. The program, as previously discussed, will not 
be able to find a converged solution in the vicinity of the snap buckling 
load . 

B. BIFURCATION BUCKLING 

The program, as currently formulated, cannot be used to solve the 
classical eigenvalue problem for the bifurcation load of an axisym- 
metrically loaded shell. However, it is possible to introduce a small 
asymmetry in the load and to determine the unsymmetric buckling behavior 
for the nearly axisymmetric load. In the limit, as the asymmetry of the 
load approaches zero, the asymmetric buckling load should approach the 
bifurcation load. This conclusion is based upon the following two mode 
analysis . 

Consider that the applied load consists of an axisymmetric part 
{, <0) } aid a very small asymmetric part It can be shown that 

due to the load and modal coupling there will be displacements in the 
modes 0, N, 2N, 4N, 6N, etc.. For this analysis only the two lowest modes 
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Snap 



Fig. 2. Typical load displacement curves for 
a shell that is sensitive to load 
imperfections . 
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are considered to be significant. According to equation (1), the 
matrix equations for these two modes are 


K 0) ] {* <0) } ■ {« (0) } - {<£>} 
K*°] - {« w } - {«ff} 


(2a) 

(2b) 


As a result of modal coupling jq^^ j- and are functions of both 

{ Z (0) } and The behavior of the axisymmetric mode is illustrated 

by curve F in Fig. 2. Note that the difference between curve E and 
curve F is due to the presence of jz in the nonlinear terms of 
equation (2a). Point P on curve F represents the maximum load for 
which a converged solution of equation (2b) can be obtained. 

The nonlinear terms in equation (2b) can be given in the form 


{«£?} - [*£?] { z(N) } 


such that 


{z ( «} - {z 


(N) 


+ Z 


(N)J 



(3) 

} } on ly . 

The solution vector 

{z<»>} 

t {z 

D 

and a nonlinear part 


7 (N) 

nL 

} 


(4a) 


where R N) } is defined by 

W 0 ] k 00 } - {* (N) } 


(4b) 


Substituting equations (3), (4a) and (4b) into equation (2b) leads to 

K N) ] {*£?} ■ Kfl K B) + 4?} <5a) 

When -^q^^j- is approximately zero, is very small. 
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If 



(5b) 


then equation (5a) becomes 

[< N) ] (Cl * [■£?] {*£?} 


Equation (5c) is an approximation to the formulation of an eigenvalue 
problem and the displacements contained in define the buckling 

mode shape. The magnitude of -jz^^ j- is negligible until the applied 
load approaches the value at point P in Fig. 2. Thus, when is 

negligible, curve F will nearly coincide with curve E and P will approach 
P' , the bifurcation point. 

In some instances, a linear solution for the axisymmetric mode is 
desirable. This is usually referred to as a linear prebuckling analysis. 
For this condition jz^j- is a linear function of the applied load, and 

, which is a function of ^Z^°^j- only> is also a linear function 
of the applied load. Therefore, equation (5c) becomes an approximation 
to a linear eigenvalue problem. Curve C in Fig. 2 illustrates a typical 
linear axisymmetric behavior with the bifurcation point denoted by point 
P' ' . Several changes had to be made to the program in order to obtain a 
linear solution in the axisymmetr ic mode and a nonlinear asymmetric 
solution. These changes are presented in Appendix B. 


thus 


k (N) 1 


C . IMPERFECTION SENSITIVITY 

An imperfection sensitive shell is one for which the buckling load 
is reduced due to the presence of asymmetric imperfections’*'. Therefore, 


1 


The sensitivity of a shell is a function of the 


loading condition. 
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curve F of Fig. 2 is a typical illustration of an imperfection 
sensitive shell since P falls below P'. For a shell that is insensi- 
tive to imperfections, P is above P f . 

In a two mode analysis, the presence of the asymmetric response will 
have some influence on the response of the axisymmetric mode since it 
appears in • «- axisymmetric mode will affect the response in 

the asymmetric mode through Thus, the sensitivity of a given 

shell and loading condition must be examined in two parts; 


(1) 

the nature 

of 

the 

effect of 

{2 « , } 

on j 

r z ( °)} 

(2) 

The nature 

of 

the 

effect of 

{z (0) } 

on j 

i* m } 


As a consequence of this interaction, the program cannot be used to 
determine imperfection sensitivity when a linear axisymmetric solution 
is employed since item (1) is eliminated. 

Application of the program to determine the sensitivity of a shell 
is limited to the extent that the asymmetric imperfections can only be 
introduced by applying an asymmetric load. Thus, the imperfections are 
not stress free as is usually assumed; nevertheless, if the asymmetric 
load is small the associated stresses will be small. 

D. SELECTION OF SHELL GEOMETRY 

In order to establish the ability of the program to determine 
bifurcation loads and the sensitivity of the shell to these loads, the 
results from the program must be compared with previously published 
analytical and experimental data. Thus, the spherical cap was selected 
as one shell geometry to investigate since numerous results, both 
theoretical and experimental, are available for a wide variety of 
loading conditions. Furthermore, previous investigations indicate 
that the spherical cap exhibits all of the traits desired; that is. 
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the symmetric behavior is a snap buckling phenomenon, bifurcation 
loads exist, for various size caps, below the symmetric snap load, and 
the shell is sensitive to imperfections over a wide range of loading 
conditions . 

Another type of spherical shell was selected for which new results 
were obtained. This was the truncated hemisphere under axial tension. 
The experimental results that are available for this shell do not agree 
with the theoretical results. Imperfection sensitivity was suspected to 
be the cause of the discrepancy. 
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IV. SPHERICAL CAP DESCRIPTION AND RESULTS 

A. GEOMETRY 

The geometry of the spherical cap is shown in Fig. 3 and can be 
specified by a nondime ns ional geometric parameter X, where 

X = 2 [3(1-v 2 )] 1/4 (H/h) 1/2 (6) 

and H is the rise of the cap, h is the thickness and v is Poisson's 

ratio. For all shells considered: 

Radii of curvature, R , R = 250 in 
5 s o 

Thickness, h = .25 in 

Modulus of elasticity, E = 30x10^ psi 
Poisson's ratio, v = .3 

B . LOAD DESCRIPTION 

The classical buckling pressure of a complete sphere is 



A uniform load q was applied to the shell from the pole out to 
r as shown in Fig. 3. 

A nondimensional load parameter X was used to describe the area 
of loading, where 

X = [ 12 (1-v 2 ) ] 1/4 (7/h)(h/R g ) 1/2 (8) 

If the approximation H = r o /2R g is substituted into equation (6), then 
equation (8) can be expressed as 

X « f- \ (9) 

o 
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For a small finite area load, approximating a point load, a 


dimensionless load parameter P* was used, where 

PR 


P* = 


and 


Eh~ 


„ -2 
P = TT r q 


( 10 ) 


(ID 


where r is equal to the distance between the pole and the first finite 
difference station. 

2 . Asymmetric 

An asymmetric load was defined in one mode as 

(N) 

q = eq (12) 

where e is a very small number. Thus, the nearly axisymmetric load 
is given by 

q T = qfl + e cos (Ne) ] (13) 


C. BOUNDARY CONDITIONS 

The first station was at the pole and was handled internally within 
the program. The final edge was clamped and the boundary conditions are 

U = V = W= § =0 (14) 

s 

where $ is the rotation of the meridian, 
s 

D . PROCEDURES 

The axisymmetric and asymmetric behavior of the spherical cap was 
determined for both uniform pressure and partial uniform loading for 
several values of X . In the majority of cases, forty stations were used 
along the meridian. The exception to this was the investigation of the 
axisymmetric snap buckling due to a uniform pressure. For 10 < X < 16, 


80 stations were used; and for X > 16, 100 stations were used. A 
stringent convergence criterion of .002 was used in all cases. 

1. Uniform Pressure 

For uniform pressure, X = X. The symmetric behavior was examined 
first in order to establish the axisymmetric snap buckling pressures over 
a wide range of X. One value of X was then selected for full examination 

of the behavior under asymmetric load. The effect of varying £ over a 

wide range was examined and the degree of sensitivity estimated. A 
value of g was selected for use in subsequent analysis of the bifurcation 
pressures. The final loads approximating the minimum bifurcation pres- 
sures were obtained for various values of X. 

2 . Partial Uniform Loading 

An investigation was made for the minimum value of X, approxi- 
mating a point load, over a wide range of X. The final loads were then 

obtained for axisymmetric snap buckling at various values of X for one 
value of X. Once the symmetric behavior had been obtained, very small 
asymmetric loads were introduced and an approximation obtained for the 
bifurcation loads as a function of X at one value of X. For small values 
of X, a linear symmetric solution was used. 

E . RESULTS 

1. Uniform Pressure 

a. Axisymmetric Snap Buckling 

The load-displacement curve at the 24 t * 1 station for X=8 is 
shown in Fig. 4. The transverse displacement at station 24 was the 
maximum displacement, and the nonlinear effects at this station were 
representative of the overall behavior. The load-displacement curve 
indicates the existence of a nearly zero slope which is indicative of 
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( 0 ) 

W /h 


Fig. 4 Axi symmetric displacement at 

station 24 versus load for a 
clamped spherical cap, \ = 8 
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snap buckling. Consequently, the final load is considered to be the 
critical pressure for axisymmetric snap buckling under a uniform 
pressure for \=8. 

Figure 5 shows the final loads obtained for several values 
of X. The results of Weinitschke [Ref. 4] and Huang [Ref c 5] are also 
presented in Fig. 5. 

b. Bifurcation Buckling 

The effect of the nearly axisymmetric load given by equation 
(13) for N=4 was obtained for 1=8. According to Refs. 4 and 5, mode four 
is the critical mode number for this case. Figure 6 shows the final load 
as a function of g. Note that as g approaches zero, the final load 
asymptotically approaches a limiting value. Figure 7 shows the dis- 
placement in mode four at station 24 versus the load for g = .0005. The 
curve strongly indicates a region of zero slope at the final load and 
convergence failure is attributed to the asymmetric behavior. Further- 
more, the nonlinear effects are definitely limited to a narrow load 

range near the final load. The ratio of nonlinear to linear displace- 

2 

ments at the final load were found to be of the order of 100 to 1, and 
the linear displacements were negligible. Thus, all of the requirements 
for a bifurcation analysis given in Section III.B. are satisfied, and 
the final load for g = .0005 is considered to be a close approximation 
to the bifurcation load. Figure 8 is a plot of the normalized, non- 
linear modal displacements versus the nondimens iona 1 , normal distance 
from the axis and is the buckling mode shape. Also shown in Fig. 8 is 
the buckling mode shape obtained by Huang for X=9 . 


The final linear displacements were computed using the initial solution 
at the first load step. 


18 




X=2[3(I-„*)] J [*]* 


Fig. 5 Critical pressures for axisymmetric snap buckling 
of clamped spherical caps. 





Fig. 7 Asymmetric modal displacement, N=4, 

at station 24 versus load, X=8, e = -0005 
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When the critical mode number is unknown, the number of 

computer runs necessary to determine the critical mode number and the 

bifurcation load can be reduced by introducing additional modes in the 

solution. For example, for \=8, an asymmetric load was introduced in 

mode one, with e = .05, and a solution was obtained in modes zero, one, 

two, three and four. Figure 9 shows the resultant asymmetric modal 

displacements at station 24 as a percentage of the total displacement 

versus the load. Note that the same softening characteristic present 

3 

in Fig. 7 exists for modes three and four . The failure of the solution 
to converge can be attributed to either mode three or four, since a zero 
slope is indicated by both. The final load is within 5 $ of the value 
obtained using a two mode solution. Thus, Fig. 9 illustrates the ability 
of a miltimode analysis to obtain a first cut approximation to the bi- 
furcation load and critical mode. 

Results were obtained for approximate bifurcation loads at 
several other values of X by selecting the critical modes as given in 
Ref. 5 and using a two mode solution with 0 = .0002. Figure 10 gives 
the final loads and critical mode numbers versus X . Also shown in Fig. 

10 are the bifurcation loads obtained by Huang, 
c. Imperfection Sensitivity 

As discussed in Section III.C., an imperfection sensitive 
shell is one for which the buckling load is reduced due to the presence 
of asymmetric imperfections. Figure 11 presents the load-displacement 
curves for the axisymmetric mode at station 24 for X=8 and several 
values of 0 . Note that as e is increased the final load is decreased. 


3 

According to Refs. 4 and 5, the bifurcation load for mode three is only 
slightly higher than that for mode four when \=8. 
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Fig. 10 Critical pressures for axisymmetric snap buckling and bifurcation 
buckling of clamped spherical caps 
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Fig. 11 Sensitivity of spherical cap, A.=8, 
to load imperfections 
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Thus, this shell appears to be imperfection sensitive for the uniform 
pressure loading condition. Fitch and Budiansky [Ref. 6] have con- 
ducted a study of the buckling and postbuckling behavior of spherical 
caps, subjected to centrally distributed pressure, based on Koiter's 
initial postbuckling theory. They found the shell X=8 to be imper- 
fection sensitive to the uniform pressure. They also predicted the 
initial slope of the bifurcation branch of the equilibrium path on a 

load-displacement plot. The predicted slope for the shell \= 8 , 

4 

subjected to uniform pressure, is given in Fig. 11 . 

2 . Partial Area Loading 

a. Axisymmetric Snap Buckling 

The final loads obtained for the axisymmetric behavior of the 
spherical cap under a small finite area load, approximating a point load 
at the pole, are shown in Fig. 12. Experimental results given by Penning 
and Thurston [Ref. 7] and Penning [Ref. 8] for axisymmetric snap buckling 
due to a small finite area load are also presented in Fig. 12. The 
experimental load-displacement curves for X > 15 show a well-defined, 
abrupt discontinuity in the displacement at the pole, and the slope of 
the curve after the snap is somewhat less than the initial slope. Below 
a value of ^. = 15, the experimental results show no axisymmetric snap"*, but 
a large majority of the load-displacement curves show significant decreases 


The load-displacement plot used by Fitch and Budiansky is based on an 
average deflection parameter. The load-displacement plot shown in 
Fig. 11 is based on the maximum displacement. The assumption is made 
that the rate of change of the average deflection parameter is essentially 
the same as the rate of change of the maximum displacement. 


"*Three atypical shells with poor radial geometry did buckle at loads some- 
what below P*=2.0, [Ref. 8]. 
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Fig. 12 Final loads for the axisymmetric behavior of clamped 

spherical caps under very small finite area loading at the pole. 



in slope in the vicinity of P*=2.0. Figure 13 shows the theoretical 
load-displacement curve for \=12 , and a reasonably accurate reproduction 
of the experimental results given in Ref. 8 for X=12.56 and a small 
finite area load. Note that although the solution failed to converge 
for P* > 2.27 snap buckling apparently is not imminent since the experi- 
mental results indicate that no snap occurs at this value of load, and 
the theoretical load -displacement curve does not show the significant 
decrease in slope that appears in Figs. 4 and 7 for the pressure loaded 
cap . 

The results given in Fig. 12 indicate that snap buckling or a 
significant decrease in the slope occurs at approximately the same value of 
load for all values of X considered. There appears to be a reasonable 
explanation for this since the experimental and theoretical results show 
the region of buckling to be in the neighborhood of the pole. For large 
values of X, the outer boundary has very little effect on the behavior at 
the pole. As X is reduced, the outer boundary moves closer to the pole 
and appears to reduce the severity of the snap, until finally, when 
X < 15, it disappears. However, the localized nonlinear effect is still 
present and shows up as a decrease in slope of the load-displacement curve. 

The final loads due to axisymmetric behavior for X=12 and 
various values of X are given in Fig. 14. As X increases, and the load 
covers more of the shell's surface, the load-displacement behavior 
exhibited the softening behavior typical in Figs. 4 and 7. 
b. Bifurcation Buckling 

The final loads for asymmetric behavior for \=12 and various 
values of X are presented in Fig. 14 with the results obtained by Fitch 
and Budiansky. The load is given by equation (13) with e=-0002. The 
critical modes were selected from Ref. 6. For X < 2.0 the final loads 
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Fig. 13 Displacement at the pole of a clamped spherical cap, 

\ = 12 , subjected to a small finite area load versus load. 
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A Asymmetric 
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a linear symmetric analysis 
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Fig. 14 Final loads for the axisymmetric and 
asymmetric behavior of a clamped 
spherical cap, X=12 , under finite 
area loading as determined by \ 





were obtained using a linear solution in the axisymmetric mode since 
the nonlinear axisymmetric solution failed to converge for loads much 
smaller than the bifurcation loads obtained by Fitch and Budiansky. 
c. Imperfection Sensitivity 

Results for the imperfection sensitivity of the shell \=12 
given in Ref. 6 show that for \ > 5 the shell is imperfection sensitive, 
but for X < 2 the load increases after bifurcation and hence, the shell 
is insensitive. This latter feature is supported by the experimental 
results given in Refs. 7 and 8 where shells were observed to transition 
into the asymmetric mode shape without a noticeable loss of load carrying 
capacity . 

It was not possible to use the program to determine the 
sensitivity for X < 2 since the solution failed to converge due to 
strong nonlinear behavior in the axisymmetric mode at loads much lower 
than the bifurcation loads . 


As discussed in Section III.C., the use of a linear solution in the 
axisymmetric mode eliminates the ability of the program to determine 
sensitivity. 


32 



V. TRUNCATED HEMISPHERE DESCRIPTION AND RESULTS 


A. GEOMETRY 


Figure 15 shows the geometry of the truncated hemisphere. The values 
of the parameters used to describe two different shells are as follows: 


Radii of curvature , R , R 

s ( 

Thickness , h 
Included angle, 0 q 
M odulus of elasticity, E 
Poisson's ratio, v 


Shell A 
120 in 
.24 in 
23.3 deg 


Shell B 
151 in 
.10 in 
23.2 deg 


30x10^ psi 
.3 


B. LOAD DESCRIPTION 

1 . Axisymmetric 

As shown in Fig. 15, an axisymmetric axial load P o /cos was 
applied to the top edge of the shell. For this loading condition, the 
reaction load at the lower boundary is P . 

2 . Asymmetric 

Two types of asymmetric loads were used; one was a small pertur- 
bation in the axial load 6P cos (Ne)/cos 0 , and the other was a small 

o o 

uniform pressure e q Q cos(Ne), where q^ is given by equation (7). Since 
the axial load was introduced as a boundary condition, and since the 
program considered the boundary conditions to be the same in every mode, 
some modifications had to be made to the program to allow a 6 different 
from 1. The appropriate statements that allow a variable 6 are presented 
in Appendix C. 
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\ P 0 (REACTION) 


Fig. 15 Geometry of the truncated hemisphere 


3b 


BOUNDARY CONDITIONS 

The boundary conditions imposed on the shell at the top edge were 


N cos 0 - Q sin 0 = P 

s o s o o 


U sin 0 - W cos 0 =0 

o o 


V = 0 


$ = 0 
s 


(15a) 

(15b) 

(15c) 

(15d) 


where N is the membrane force per unit length, and Q is the transverse 
s s 

force per unit length. 

At the reaction end the shell was clamped, and the boundary conditions 
are given by equation (14). 


D . PROCEDURES 

A forty point finite difference net was used in all cases, and the 
convergence criterion was .01. The nonlinear axisymmetric behavior for 
the two shells was first obtained. Both types of asymmetric loads were 
then introduced, and the final loads, mode shapes and load-displacement 
curves were obtained for a linear and a nonlinear solution in the axisym- 
metric mode. The sensitivity of shell B was examined by varing both 6 
and e • 


E . RESULTS 

1. Axisymmetric Behavior 

The linear and nonlinear parts of the radial displacements at the 
final load are shown in Fig. 16 for shell A. The final displacement state 
of shell B was similar. Figure 17 depicts the load-displacement behavior 
for a station near the top of the shell where maximum stiffening occurred. 
Since Fig. 16 shows that there are no softening nonlinearities, failure 
of the solution to converge is attributed to the stiffening nonlinear 
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W (o) /h 


Fig. l6 Axisymmetric , normal displacements 

of Shell A under axial tension 

P /Eh=l. 325x10-3 
o 
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behavior. The final nondimensional load P /Eh was 1.325x10 for shell 

o 

-3 

A and .600x10 for shell B. 

2 . Bifurcation Buckling 
a. Shell A 

The theoretical results obtained for the bifurcation loads 
versus mode number and the experimental buckling load given by Yao 
[Ref. 9] are presented in Fig. 18. Yao's theoretical results are based 
on a linear, membrane prebuckling assumption. Also shown in Fig. 18 are 
the present results obtained for the final loads using a linear symmetric 
analysis with an asymmetric load in the appropriate mode and 6=.0001. 
Figure 19 presents the buckling mode shape for N=40 and 6=.0001, and 
the mode shape assumed by Yao. The linear portion of the displacements 
was negligible at all stations along the meridian, and the requirements 
for a bifurcation analysis are satisfied. Similar results were obtained 
for N > 30. 

For N < 30, the results for the final loads differ signifi- 
cantly from the bifurcation loads predicted by Yao. Figure 20 shows the 
modal displacements along the meridian at the final load for N=15 and 
6=.0001. Examination of Fig. 20 reveals that the linear part of the 
solution is larger than the nonlinear part, and that there is stiffening 
in the region near the applied load. Hence, the solution does not satisfy 
the conditions for a bifurcation analysis, and therefore, the final loads 
obtained for the lower mode numbers do not represent a critical load 
condition . 

The values of the final loads were also obtained using a 
nonlinear symmetric solution. Figure 21 shows the nonlinear and linear 
parts of the modal displacements at the final load for N=40 and 6=.0001. 
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Linear Symmetric 
Analysis 


Convergence Failure Due To Symmetric Nonlinear 
Behavior 


Experimental Bucking Load [Ref. 9] 
(Mode Number Not Given) 

I 


MODE NUMBER 


Final loads for asymmetric behavior of 
Shell A subjected to axial tension with 
small asymmetric load, 6 = .0001 






Fig. 20 Meridional mode shape, N=15, of 

Shell A subjected to axial tension 
with small asymmetric load, 6=.0001 


4l 
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Fig. 21 Development of buckling mode shape, N=40, 
of Shell A subjected to an axial tension 
load in the vicinity of the critical load 
with small asymmetric load, 6=.0001 


1+2 




I 

Note that a mode shape is emerging, but the conditions for bifurcation 
analysis are only satisfied over that part of the shell where the linear 
displacements are negligible. Furthermore, convergence failure could 
not be definitely associated with either the asymmetric or the axisym- 
metric behavior since the final load for axisymmetric nonlinear behavior 
and the bifurcation load are approximately equal, 
b. Shell B 

The axisymmetric, nonlinear behavior of shell B did not 
experience convergence failure until the load was much larger than the 
minimum bifurcation load predicted by Yao. As a consequence, results 
were obtained for the final loads using a linear and nonlinear symmetric 
analysis with an asymmetric load in the appropriate mode and 6=.0001. 

The results are presented in Fig. 22 along with the theoretical bifur- 
cation loads and experimental buckling load presented in Ref. 9. 

Figure 23 shows the nonlinear part of the modal radial 
displacements at the final load for N=70, e = -001, and a nonlinear 
symmetric analysis. The linear part of the modal displacements were 
negligible at all stations and the mode shape is essentially the same 
as the mode shape for shell A shown in Fig. 19. Although the solution 
was obtained using a nonlinear symmetric analysis, failure of the 
solution to converge could definitely be associated with the asymmetric 
behavior . 

For N < 35 final loads were obtained below the bifurcation 
loads predicted by Yao. As was the case for shell A, the solutions at 
these lower mode numbers did not satisfy the conditions for a bifurcation 
analys is . 
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Nonlinear Symmetric Analysis 
Linear Symmetric Analysis 


Experimental Bucking Load [Ref. 9 ] 
(Mode Number Not Given) 








3 . Imperfection Sensitivity 


The sensitivity of shell A could not be examined since a linear 
axisymmetric solution was required for the bifurcation analysis. 

The sensitivity of shell B for N=70 was examined using an asym- 
metric axial load 6 P q Cos (N e)/cos 0^ and an asymmetric pressure load 
eq o Cos(Ne). The variation of final load with respect to 6 and e is 
shown in Fig. 24. For the asymmetric axial load, the final load was 
constant as the magnitude of § was increased. However, examination of 
Fig. 21 reveals that the asymmetric axial load causes linear displace- 
ments only in the region near the applied load. Outside this region the 
linear displacements are effectively zero^ . Therefore, the use of an 
asymmetric axial load does not introduce linear displacements of any 
consequence away from the region near the applied load. On the other 
hand, the use of an asymmetric pressure load introduces linear displace- 
ments over the entire region of the shell and is considered more appli- 
cable for examining sensitivity. 

Figure 25 shows the effect of the asymmetric pressure load upon 
the final load and axisymmetric displacement at station twenty. The 
behavior at all stations in the central region of the shell is exemplified 
by the behavior at station twenty. For e < .01, there was a small decrease 
in the final load and a small increase in the final axisymmetric displace- 
ment. However, for e=.10, there was a significant reduction in the final 
load, and the final axisymmetric displacement is approximately the same 
as that associated with e=.01. This type of behavior is similar to that 
obtained for the spherical cap under uniform pressure. 


Since the net axial load in the asymmetric mode is zero, the linear modal 
stresses become smeared as the distance from the load application is 
increased. Thus, there is little or no asymmetric stress present at the 
reaction end . 
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Fig. 24 Variation of final load with 6 and e for Shell B 



Fig. 25 Sensitivity of Shell B to load imperfections 
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Hutchinson [Ref. 10], in an analytical study of the initial 
postbuckling of toroidal shells, suggested that the imperfection 
sensitivity of Yao's truncated hemispheres should agree qualitatively 
with the degree of sensitivity that he obtained for the simply sup- 
ported hemispherical segments subjected to axial tension. Shells A and 
B fall within the imperfection sensitive range predicted by Hutchinson. 

Q 

A plot of the initial slope of the generalized load-deflection curve , 
as given in Ref. 10, for the perfect toroidal shell corresponding to 
shell B is shown in Fig. 25. 


g 

Hutchinson's load-displacement plot is based on axial elongation. An 
appropriate conversion was applied to his slope to account for differ- 
ences in nondimensionalization. 
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VI. CONCLUSIONS 


A. FINAL LOADS AS BUCKLING LOADS 

t 

As a consequence of the method of solution used in the computer pro- 
gram, the solution may fail to converge for both softening and stiffening 
nonlinearities. Thus, a careful examination of the solution should be 
made before proclaiming a final load to be a buckling load. For axisym- 
metric snap buckling, the load-displacement curves should indicate that 
a region of nearly zero slope has been reached. For bifurcation buckling 
loads, the above should hold true for the modal displacements, and the 
linear part of the modal displacements must be negligible compared to 
both the nonlinear part of the modal displacements and the axisymmetric 
displacements . 

There are instances when it's not completely clear why the solution 
fails to converge. For example, consider the results obtained for very 
small finite area loads on a clamped spherical cap. For X > 15, the 
experimental results indicate that an axisymmetric snap occurs at an 
approximate value of load for which the computer solution failed to 
converge. However, for \ < 15 the solution experienced convergence 
problems at loads where the slope changed significantly, but no snap 
buckling occurred. Archer [Ref. 11] obtained similar theoretical results 
for simply supported caps using the pseudo load approach. Mescall [Ref. 
12] notes that Archer's final loads were loads for which the numerical 
procedure could not converge and were not due to the existence of a 
maximum on the load-displacement curve. Throughout this study an attempt 
has been made to clarify this point. In those cases where the results 
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did not indicate the existence of a maximum on the load-displacement 
curve , the final loads were not classified as buckling loads. 

B. AXISYMMETRIC SNAP BUCKLING 

1 . Spherical Cap 

The results obtained for the axisymmetric snap buckling loads of 
the spherical cap under uniform pressure compare favorably with previ- 
ously published theoretical results. Since the load-displacement curve 
shows the approaching region of zero slope associated with the existence 
of a maximum and snap buckling, and since convergence failure can definitely 
be attributed to this softening nonlinear behavior, the classification of 
the final loads as snap buckling loads is considered justified. 

The case of small area loading, approaching a point load, gave 
results for high rise caps that compare favorably with experimental 
results. The load-displacement behavior indicated that snap buckling 
was iminent . At smaller values of X, the behavior is not characterized 
by indications that a maximum exists on the load-displacement curve, and 
failure of the solution to converge is not associated with snap buckling, 
but with moderate changes in stiffness of the shell. 

As the area of loading is increased, the behavior is again 
characterized by the type of nonlinear behavior indicative of snap 
buckling and the final loads obtained are considered to be snap buckling 
loads . 

2 . Truncated Hemisphere 

Final loads were obtained for the axisymmetric behavior of the 
truncated hemisphere under axial tension. However, these loads do not 
represent critical loads, but are the last value of the load for which 
a solution could be obtained. Failure of the solution to converge is 
attributed to stiffening nonlinearities. 
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C. BIFURCATION BUCKLING 


The results of the investigation of asymmetric behavior under nearly 
axisymmetric load indicate that the program can be used to predict 
bifurcation buckling loads. In fact, the results obtained for the 
bifurcation loads using the asymptotic procedure are better than was 
anticipated, since an extremely small disturbance in the mode is all 
that is required to trigger nonlinear modal displacements that are much 
larger than the linear modal displacements as the load approaches a 
critical value. 

1 . Spherical Cap 

The bifurcation loads obtained for the spherical cap for all 
loading conditions are in very good agreement with previously published 
theoretical results. For an asymmetric load with a maximum magnitude of 
.01$ to .05$ of the axisymmetric load, the linear asymmetric displacements 
are negligible compared to the axisymmetric displacements and the nonlineal 
asymmetric displacements are negligible until the applied load reaches a 
value within 95$ of the final load. The buckling mode shapes compare 
favorably with previous results. 

2 . Truncated Hemisphere 

The bifurcation loads obtained for the truncated hemisphere are 
in good agreement with the previously published theoretical results, and 
the critical mode numbers are in excellent agreement. Several final loads 
were obtained for low mode numbers that do not agree with the previous 
results, but these loads are easily eliminated as bifurcation loads. 

D. IMPERFECTION SENSITIVITY 

The sensitivity of a given shell and loading condition to imperfection.': 
is a result of the nature of the modal coupling between the axisymmetric 
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mode and the asymmetric modes. When the failure of the solution to 
converge can definitely be attributed to softening nonlinearities in 
the asymmetric mode, the computer program can be used to determine the 
coupling effects and the sensitivity can be estimated. 

1. Spherical Cap Under Uniform Pressure 

For the pressure loaded cap with X=8, reductions in the final 
load with increasing magnitude of imperfections are accompanied by 
reductions in the final value of the axisymmetric displacements. Further- 
more, the failure of the solution to converge is due to softening non- 
linearities in the asymmetric mode. Therefore, this pressure loaded cap 
shows a strong sensitivity to imperfections. This conclusion is in 
agreement with the previous theoretical analysis given in Ref. 6. 

2 . Truncated Hemisphere Under Axial Tension 

The axisymmetric behavior of the truncated hemisphere under axial 
tension shows stiffening over the entire region of the shell, and the 
solution encounters convergence problems due to this stiffening. Never- 
theless, the final load for shell B, subjected to a very small asymmetric 
load in the critical mode, is a bifurcation load, and a sensitivity analysis 
of this shell is possible. For increasing values of asymmetric axial load 
there is no change in the final load. This is due to the fact that the 
imperfections in the axial load generate displacements and stresses only 
in the upper region of the shell. In this region, the asymmetric dis- 
placements have negligible effect upon the axisymmetric displacements, 
and the shell appears to be imperfection insensitive. On the other hand, 
an asymmetric pressure load causes displacements over the entire shell. 

For very small values of e, less than .01, the asymmetric mode does have 
some effect on the axisymmetric mode; the final load is slightly reduced 
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but the final displacements in the axisymmetric mode are increased. 

Thus, the shell is sensitive to small imperfections, but the behavior 
differs from that of the spherical cap. For larger imperfections the 
asymmetric mode Ijas a pronounced effect upon the axisymmetric behavior 
in the central region of the shell, and there is a significant reduction 
in the value of the final load with virtually no change in the final 
value of the axisymmetric displacements. This type of behavior is 
similar to that of the spherical cap. 

Hutchinson [Ref. 10] suggested that the imperfection sensitivity 
of Yao's shells should agree qualitatively with the degree of imperfection 
sensitivity he obtained for simply supported hemispherical segments. 
Reference 13 contains a summary of the results of Hutchinson's study and 
a discussion on the postbuckling coefficient used to determine sensitivity. 
The imperfection sensitivity of circular cylindrical shells under axial 
compression is used as a calibration for the significance of the post- 
buckling coefficient. Based on Hutchinson's prediction, shell B is not 
as sensitive to imperfections as the cylinder or the spherical shell, 
but it does fall in a region where a significant reduction in the buckling 
load may occur . 

The present results for small e are in agreement with Hutchinson's 
prediction. However, for larger e the results indicate an apparent change 
in the nature of the sensitivity and the behavior resembles that of the 
spherical cap. Since the experimental buckling loads of shells A and B 
are about one half of the theoretical bifurcation loads, which is common 
for pressure loaded spherical caps, this analysis appears to offer a 
reasonable explanation for the large discrepancy between the bifurcation 
loads and the experimental buckling loads. 
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APPENDIX A 


USER PREPARED SUBROUTINES 

The following are user prepared subroutines necessary to describe 
the stiffness, geometry and loads for the various shells examined in 
this study. Certain input parameters were required before calculation 
of any quantities in the subroutines. These values were input into the 
main part of the program. 

KMAX - total number of meridional stations 

MNMAX - the number of terms used in the expansion of the 
applied load 

NU - Poisson's ratio (v) 

EIAST - reference modulus of elasticity (E ) 

TKN - reference thickness (h ) 

CHAR - characteristic length (a) 

SIGO - reference stress (a ) 

o 


A. STIFFNESS 


The stiffness quantities necessary in the program are calculated at 
each meridional station (K) by subroutine BDB. The nondimens ional form 
of the Fortran variables for a homogeneous, isotropic, constant thickness 
shell are as follows: 


I s l 

E h (1-v 2 ) (1-v 2 ) 

o o 


(E = E , h = h ) 
o o 


/e C 2 dC 

E h 3 (1-v 2 ) 
o o 


B_ 

12 


DB = f = 0 
ds 


DD 


dD 

ds 


0 
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Since all the shells examined were assumed to be homogeneous, 
isotropic and constant thickness, only one subroutine was required. 


SUBROUTINE BDB (K,B ,DB ,D,DD) 

REAL NU 
COMMON 

1/BL32/TKN,EIAST,CHAR,SIG0/BL15 

2/NU,Ul (10) ,V1 (10) ,W1 (10) , V2 (10) ,U2 (10) ,W2 (10) ,U3 (10) , 
3V3 (10) ,W3 (10) 

B=1 . / (1 . -NU**2) 

D=B/12 . 

DB=0 . 

DD=0 . 

RETURN 

END 


B . GEOMETRY 

The geometric quantities necessary in describing the shells are 
defined in subroutine GEOM. The following are the nondimensional 
Fortran variables required at each meridional station (K) . 


DEL = 


(KMAX-l)a 

R(K > ' [f] K 
GAM(K) = [££!•] K 

0MI(K > ■ [r] K 

o"* 1 ® - [r] k 

s 

DEOHX (K) - ^ [|-] 


K 


1 . Spherical Cap 

Prior to calling subroutine GEOM, certain quantities were input 
or calculated in MAIN and are as follows: 
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FCTZ = X 

CHAR = R = R 

HTR = H S ® 

FCTL = S 

RL = r 


(input) 
(input) 
(calculated) 
(calculated) 
(ca lculated) 


The necessary statements in MAIN required for calculating the desired 


quantities are as follows : 

HTR=TKN* (FCTZ **2 ) /SQRT (5 8 . * ( 1 . -NU**2 ) ) 
RL=SQRT ( (2 . *CHAR*HTR) - (HTR**2 ) ) 
FCTL=CHAR*ARS IN (RL/CHAR) 


Subroutine GEOM is as follows : 


SUBROUTINE GEOM 
COMMON 

1/IBL4/KMAX,KL 

0/BL8/R(200) ,GAM(200) ,OMT(200) 
3/BLll/OMXI(200) ,PHEE,T0,T2 
4/BL17/DEL 
4/BL20/DEOMX (200) 

6/BL32 /TKN , ELAST , CHAR , S IGO 
RK*=KMAX- 1 

DEL= (FCTL/RK) /CHAR 
DO 4 K-2 , KMAX 
RK=K 

thet=(rk-i.)*del 

R (K) =S IN (THET ) 

GAM(K) =COS (THET)/R(K) 
OMT(K)=l. 

OMXI (K)=l . 

4 DEOMX (K)=0 . 

R(1)=0. 

GAM(1)=0 . 

OMT ( 1 ) = 1 . 

OMXI (1)=1 . 

DEOMX(1)=0. 

RETURN 

END 


2 . Truncated Hemisphere 

The input quantities required in MAIN are as follows: 

PHIO = 0 (input in radians) 

CHAR = R° = R (input) 

S 0 

Subroutine GEOM is as follows : 


57 



SUBROUTINE GEOM 


(COMMON STATEMENTS) 

RK=KMAX-1 
DEL=PHIO/RK 
DO 13 K=1 ,KMAX 
RK=K-1 

BETA=PHIO- (RK*DEL) 
R(K)=1 . - ( 1 . -COS (BETA) ) 
GAM(K) =SIN(BETA)/R(K) 
OMT(K)=COS (BETA)/R(K) 
OMXI(K)=l. 

13 DEOMX(K)=0. 

RETURN 

END 


C. LOAD 

The axial loading case is covered in Appendix C and the pressure 
loading case is described in the following. 

The reference pressure loads for each coefficient used in the cosine 
expansion of the load, at each meridional station, were calculated or 
defined in subroutine PIOAD. The following are the nondimens ional Fortran 
variables required 

NN(I) = n 

PR(i) - §V 

o o 

For the shells examined, the following reference values were used and 

were input or calculated in MAIN. 

MODI - the first asymmetric mode number used in load 

description (input) 

CHAR = a = R = R (input) 

TKN = h = h S ® (input) 

SIGO = ct°= q (calculated) 

o o 

The statement required to calculate SIGO is 

S IGO= (2 . *ELAST*TKN) / (CHAR*SQRT (3 . * ( 1 . -NU**2 ) ) ) 
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1 . Spherical Cap 

The load description for the spherical cap, in all cases, 
required a symmetric pressure load and, for the majority of runs, one 
or more modal pressure loads. The symmetric mode (n=0) was always 
NN(1) and PR(1) was equal to -1. corresponding to an external pressure. 

A Fortran variable, FCTR, was an input parameter determining the number 
of stations loaded and was equal to KMAX for the uniform load. The 
magnitude of the asymmetric load was governed by the input variable EPSL. 
The subroutine PLOAD for the spherical cap is as follows : 
SUBROUTINE PLOAD (K) 

COMMON/BL32 /TKN , ELAST , CHAR, S XG0/IBL2/NN (10) ,MNINIT 
8/BL6/Z (4 ,220) ,S0E , OSE ,ALQAD/IBL1/MNMAX 
5/BL3/PR(10) , PX (10) , PT (10)/IBL4)KMAX,KL 
6/BL8/R(200) , GAM (2 00) , OMT (200)/BIL2/EPSL,MOD1 , LINA 
M=M0Dl-2 
DO 10 1=1 ,MNMAX 

10 NN(I)=M4-I 
NN(1)=0 

IF(K.GT.FCTR) GO TO 39 
DO 11 1=1 ,MNMAX 

11 PR(I)=EPSL 
PR(1)=- 1 . 

GO TO 40 

39 DO 12 1=1 ,MNMAX 

12 PR(I)=0 . 

40 RETURN 
END 

2 . Truncated Hemisphere 

For the truncated hemisphere the only pressure load used was an 
asymmetric load in one mode, and PR(1) was always equal to zero. The 
Fortran parameter EPSL was used to govern the magnitude of the asymmetric 
axial load and the magnitude of the asymmetric pressure load was specified 


within subroutine PLOAD. If an asymmetric axial load was used then PR(2) 
was set equal to zero. Subroutine PLOAD for the truncated hemisphere is 
as follows : 
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SUBROUTINE PLOAD(K) 

(COMMON STATEMENTS) 

NN(1)=0 
NN(2)=MODl 
PR ( 1 ) =0 . 

PR(2)=- . 001 

RETURN 

END 


6o 



APPENDIX B 


ADDITIONS TO THE PROGRAM FOR LINEAR SYMMETRIC SOLUTION 
The following is a listing of the statements and their location 
within the various sections of the program that allowed the nonlinear 
terms in the symmetric mode to be by-passed and the linear solution 
obtained . 


f 


MAIN 


(COMMON BLOCKS) 


C0MM0N/BIL2/EPSL ,M0D1 , LINA 

t READ(5,102)NU,SIGO,EIAST,TKN,CHAR 

C SET LINA=1 FOR LINEAR PREBUCKLING 
LINA-1 


^ READ (5, 102) DELOAD , EPS 

t 10 WRITE (6,201) NO 

IF (LINA. EQ. 1)WRITE (6 , 118) 

118 FORMAT (48X, 'LINEAR PREBUCKLING'//) 


WRITE (6 ,100) TITLE 

v 


f 

t 


SUBROUTINE XANDZ 


(COMMON BLOCKS) 

COMMON/BIL2 /EPS L, MODI , LINA 


031 

032 

038 

039 

051 

052 

531 

564 


6l 


565 


I 

| DO 8 M=1 ,MNMAX 622 

CALL TLOAD(l) 623 

IF(ik.EQ. 1. AND. LINA. EQ.l) GO TO 7 


FFS(l,M)=-TT(M)*ALQAD+OSE*(BXl(M) ETC. . . 

FFS (2 ,M)=QSE*(B1 *D1 *BST1(M ETC... 625 

FFS (3 ,M)=IAM2*GAMl*Dl*MT(M)*ALOAD- . . . ETC ■ , ■ 


GO TO 8 
7 FFS(1,M)=0. 
FFS (2 ,M)=0 . 
FFS (3,M)=0. 


8 FFS (4 ,M)=0 . 

V 

I DO 14 M=1,MNMAX 

CALL TLOAD(KMAX) 

IF(M. EQ.l. AND. LINA. EQ.l) TO TO 11 


FLS (l)=-TT(M)*ALOAD+OSE*(BX3 (M) . . 

. . .ETC. . . 

651 

FLS (2 )= OS E * (B L*D 1 *BXT3 (M)+EX3 . . . 

. . . ETC . . . 

652 

FLS (3)=IAM2*GAML*Dl*MT(M)*ALOAD. . . 

. . . ETC . . . 

653 


GO TO 13 
11 FLS (1)=0. 
FLS (2)=0 . 
FLS (3)=0 . 


13 IK=KL+KMAX* (M- 1) 

\l 


654 


SUBROUTINE FORCE (K) 

(COMMON BLOCKS 


960 


987 


CQMM0N/BIL2/EPSL, MODI, LINA 


A EX2T=EX2(M) 
ET2T=ET2 (M) 


988 

066 

067 


7 IF(LINA..EQ. l.AND.M.EQ. 1) TO TO 31 
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068 


GEE ( 1 ) =GEE ( 1 ) *OSE* (BS* (DBX+DBE+GA* ETC . . . 

2 (EXT2T+ETT2T) ) ) *TDEL 

31 IF(K.GT.l) GO TO 10 
DO 20 1=1,4 

SUBROUTINE OUTPUT(IMODE) 

(COMMON BLOCKS) 

C0MM0N/BII2/EPSL, MODI , LINA 

f 

f TTS =TT (MN ) *ALOAD 

IF(MN.EQ.l. AND. LINA. EQ.l) GO TO 13 

| EX=(U3(MN)-U1(MN))*TDLI +OX*W2 (MN) + ETC .77 

^ 1 + ENL*SOE*BXT3 (MN)) 

GO TO 14 

13 EX= (U3 (MN) -U 1 (MN) ) *TDLI+OX*W2 (MN) 

ET=ENR*V2 (MN)+GA*U2 (MN)+OT*W2 (MN) 

EXT= . 5*( (V3 (MN) -V 1(MN))*TDLI-ENR*U2 (MN) -GA*V2 (MN) ) 

14 KT=ENR*PHIT(MN)+GA*PHIX(MN) 

KXT= . 5 * (ENR* ( - PHIX (MN ) -GA*W2 (MN) + (W3 ETC... 

V 


076 

078 

445 

473 

474 

521 

522 

525 

526 

527 


CALL TLOAD(K) 

TTS=TT (MN ) *ALOAD 

IF(MN.EQ. 1. AND. LINA. EQ.l) GO TO 15 

EX= (U3 (MN) -U 1 (MN) *TDLI +OX*W2 (MN) ETC... 

1 + ENL*SOE*BXT3 (MN) ) 


851 

852 


853 


856 


GO TO 16 

15 EX=(U3 (MN) -U1 (MN) )*TDLI+OX*W2 (MN) 

ET=ENR*V2 (MN)+GA*U2 (MN)+OT*W2 (MN) 

EXT=. 5*( (V3 (MN)-Vl (MN))*TDLI - ENR*U2 (MN) -GA*V2 (MN) ) 

16 KT=ENR*PHIT(MN)+GA*PHIX(MN) 857 
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APPENDIX C 


ADDITIONS TO THE PROGRAM FOR ASYMMETRIC AXIAL LOAD 
In order to adjust the magnitude of the asymmetric axial load, 
additional logic was incorporated in subroutines XANDZ and FORCE. 

The axial loads are entered as boundary conditions and are described 
by the Fortran variables ELIS and ELLS. A Fortran variable EPSL was 
used to determine the magnitude of the asymmetric load and is contained 
in the COMMON block BIL2, which was introduced in Appendix B. The sym- 
metric part of the load was always taken as the first mode, described 
by NN(1)=0. The necessary statements required are as follows: 


SUBROUTINE XANDZ 

DO 14 1=1,4 
SUMZ=0 . 

DO 15 J=1 ,4 

IF (M . EQ . 2 . AND . I . EQ . 1 ) ELLS(J)=ELLS (J)*EPSL 

15“ SUMZ=SUMZ+ZFIM(I,J,M)* ETC. . . ' 

SUBROUTINE FORCE (K) 

DO 21 J=1 ,4 

IF (M.EQ.2 .AND.I.EQ. 1) ELIS (J)=EL1S (J)*EPSL 
21 SUMX=SUMX+DL (I~ J , M) *E LIS ETC . . . 


531 

657 

658 

659 


660 

960 

081 


082 
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